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Abstract 

A detailed comparison of ten low-Reynolds-number k — e models is carried 
out. The flow solver, based on an implicit approximate factorization method, is 
designed for incompressible, steady two-dimensional flows. The conservation of 
mass is enforced by the artificial compressibility approach and the computational 
domain is discretized using centered finite differences. The turbulence model 
predictions of the flow past a hill are compared with experiments at Re — 1.4 • 
10 6 . The effects of the grid spacing together with the numerical efficiency of 
the various formulations are investigated. The results show that the models 
provide a satisfactory prediction of the flow field in the presence of a favourable 
pressure gradient, while the accuracy rapidly deteriorates when a strong adverse 
pressure gradient is encountered. A newly proposed model form that does not 
explicitly depend on the wall distance seems promising for application to complex 
geometries. 


1 Introduction 

The direct numerical simulation (DNS) of turbulent flows[13] is restricted to simple 
configurations and low Reynolds numbers (Re) because of the speed and memory lim- 
itations of supercomputers. A valid alternative, given by the Reynolds averaging tech- 
nique of the Navier- Stokes equations with a model of Reynolds-stress, allows solving 
engineering problems with a manageable computer effort. Following the Boussinesq 
assumption , the Reynolds stresses are proportional to the mean velocity strain: 


- UiUj = v t 


(dU 3 dU { \ 

\ dx'i + dxj ) 
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in which k = ~ Ya= i u i u i 1S the turbulent kinetic energy, U{ the fluctuating velocity, U{ 

the mean velocity, v t is the so called turbulent viscosity , and the overbar stands for the 
average operator. The problem now is how to determine v t . With the local velocity 
and length scale of turbulence, v t can be modeled by the Kolmogorov expression: 

k 2 

v t = c »~ ( 2 ) 

in which e = vvjJUiJ is the dissipation rate of turbulent kinetic energy. The quantities 

k and e will be determined by corresponding transport equations. 

This approach has been applied to a very wide class of two dimensional and three 
dimensional compressible and incompressible flows [20]. Most of the applications use 
the high Re form of the turbulence models in which the effect of molecular viscosity 
is neglected and the presence of solid walls is accounted for by means of so called 
wall functions[ 11]. The model equations have been tuned using a set of experimental 
channel data obtained in the absence of strong pressure gradient. 

As pointed out by Bradshaw[3], the effect of streamline curvature, induced by a 
pressure gradient not aligned with the main flow direction, is not properly accounted 
for in these models. Problems are even more evident in the case of separation since 
the wall function approach is strictly valid only for attached shear layers provided 
that the turbulence is in local equilibrium in the wall proximity. Several authors, 
among them Rodi and Scheuerer [22], showed that standard two equation models 
fail to reproduce the correct flow pattern in the case of adverse pressure gradients. 
Although the limitations of these models are well known, their better performance 
may be obtained by proper modeling of the viscous and buffer layers, which are usually 
bridged by wall functions. 

With this in mind new formulations of the two-equation models, called low Re 
forms ( LR ) have been introduced in the past several decades. These models ac- 
count for the molecular viscosity which is not negligible in wall regions, the normal 
velocity fluctuation damping induced by a solid boundary, and the presence of a 
anisotropic contribution to the dissipation rate of turbulence that dominates in the 
viscous layer[19, 24]. In the LR models the turbulence quantities are solved through- 
out the viscous and buffer layers thereby accounting for convection and diffusion in 
the molecular viscosity dominated region. The LR forms have been further tuned to 
lit the experimental data in the viscous and buffer layers, but, because of the difficul- 
ties in performing measurements in the wall proximity, the models have been often 
made to reproduce measurements with a very wide uncertainty range. Thanks to the 
recent availability of DNS results, although limited to very simple flows and quite 
low Reynolds numbers, it has been possible to formulate new LR models that are con- 
sistent with these numerical data sets. In a previous study by the authors[17], a new 
model has been compared with several other LR forms. The preliminarj r investigation 
was carried out by comparing the model predictions with DNS of a fully developed 
channel flow. In the present paper the attention will be focused on LR forms applied 
to a more complex geometry in the presence of streamline curvature and separation 
to highlight the differences in the various models. From a computational point of 
view it is important to observe that this investigation is carried out using a single 
solver able to cope with all the models thereby ensuring the same degree of accuracy 
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in the solution of the various model equations. Together with the turbulence models, 
the implicit solver will be briefly described in the following. 


2 Selection of the k — e LR models 


The investigation carried out in ref. [17] includes ten models which are: Chien[4], (ch), 
Jones and Launder[9], (j/), Nagano and Hishida[18], (nh), Coakley[6], (co), Speziale 
et al.[25], (sp), Kim[8], (fci), Rodi[21], (ro), Lam and Bremhorst[10], (/&), Shih[23], 
(s/i), Michelassi and Shih[17], (ms). The letter codes in parenthesis will be used to 
refer to each model. All these models are further tested in this study. 

It is possible to express the LR models in a unified form using the following 
nondimensional variables (the overbar stands for a nondimensional quantity): 


k = 


k 

U2 





xl = 


Xi 

L 


in which U is a typical velocity, L is a typical length, and v is the fluid molecular 
viscosity. Accordingly, the flow Reynolds number is defined as: 


Re — 


L U 

v 


Hereafter the variables will appear only in their nondimensional form so that the 
overbar will be dropped for simplicity. The turbulence model transport equations 
may be written as follows: 


It + iSr = i ^ g) + p - £ + -° + n 
m + = A (4 ( l + £) it) + «iAf ^ - «hi +E 


( 3 ) 


in which 

is the production rate, and 

jfe 2 

Vt ~ Re y (4) 

is the Kolmogorov expression for the txirbulent viscosity modified by including the 
function The extra terms D, i?, II represent the corrections to the standard 
high Re formulation. They allow balancing the transport equations (3) in the wall 
region[19]. Their detailed forms are listed in Table 2, and here, 


e = e-D (5) 

In the q - lu model ( q = y/k, u; = and the k — r model (r = ~), the transport 
equations of uj and r are introduced instead of the e equation. 

The ro two-equation and ki four-equation models are two-layer models in which 
an inner layer is defined close to the wall (t/ + < 100) where, while using a transport 
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model 


A 

/2 

ch 

1 — eaip(-.0115t/ + ) 

1.0 

l-0.22e*p(-f|) 

ji 

exp (^i) 

1.0 

1 — 0.3exp(—Rf) 

nh 

(1 - exp(-^j)) 2 

1.0 

1 - 0.3ea '<p{-R 2 t ) 

CO 

1 — exp(— .0065R y ) 

1.0 

1.0 

sp 

< 1 + M- ,anA (£» 

1.0 

t 

1-0.22 exp(-§) 

ki 

1 — exp(-.025</R t - 10~ 5 R 2 ) 

1.0 

1.0 

TO 

1 - €££(-.0198.^) 

1.0 

1.0 

lb 

(1 - ea:p(— .0165-Rfe)) 2 • (1 + ^~) 

(1+f) 2 

1 

Hi 

T 

to 

sh 

1 -exp(-ZLi<^iy +i ) 

1.0 

1-0.22 exp(-§) 

ms 

_t 

-| _ exp( — .0004- exp(l. 2 Ri,)) 
exp( — .0004) 

1.0 

1-0.22 exp(-§) 


Table 1: Damping functions in the LR forms 


equation to compute the turbulent kinetic energy, the dissipation rate is computed 
by an algebraic expression which is a function of the wall distance y and turbulent 
kinetic energy k. This choice allows removing most of the numerical problems usually 
connected with the LR forms in the wall region. However, the two-layer models have 
to properly deal with the problem of matching the inner layer with the outer layer 
in which different sets of equations hold, ki and ro use similar criteria to pass from 
the inner layer to the outer layer equations; the interface is placed at v t jv ss 35 - 36. 
Although k is solved by using the same transport equation in both the layers, the 
dissipation rate e comes from an algebraic expression within the buffer layer and from 
a transport equation in the outer layer. 


2.1 The damping functions 

The damping functions selected in the various formulations may be found in table 1, 
whereas table 2 gives all the selected forms for the extra terms D , E, n. Most of the 
models have a set of damping functions closely related to the Van Driest expression 
for the mixing length, L: 

L = y( l-e-^/26) 

where t/ + is 

= u T y Re (6) 

The turbulent Reynolds numbers based on y, R y = y/k y Re, or R t ~ are also 
used to mimic effects of the wall. In nearly all the formulations under investigation 
the decay of the dissipation rate is governed by the exponential function fi proposed 
by Hanjalic and Launder [7], according to which c 2 has a finite value at the wall. The 
sp formulation needs an additional function to further damp the value of ci in the 
wall region (see ref.[25]). The empirical constant ci is left unchanged all the way to 
the wall in all the formulations with the exception of lb where, since no extra terms 
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model 

n 

D 

E 

ch 

0 

1 w 
1 

R e ^ ex P( .5y+) 

ji 

0 

2 (dVky 

Re \ dy ) 

2 i/ t (a 2 u\ 2 
Re 2 Ul/ 2 ) 

nh 
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2 (aVk) 2 

Re \ dy ) 

-j'l -fi - f ) {?v) 2 

CO 

0 

0 

- 

sp 

0 


, j 

ki 

0 

- 

! - 1 
j 

TO 

0 

- 


ib 

0 

j o 

0 

sh 


0.05 vt jL 1 

: 0 

(dlu) 2 

RP \ dy* ) 


.U(l-exp(-y+))a k j 

j ms 


[ p -r "kj) 
n >3 i ,j 

| 

; 0 

i 

Re* \8y* ) 


Table 2: Extra terms in the LR forms 


are added in both k and e, the production of dissipation rate must be damped in 
order to balance the e equation at the wall. 

The choice of the damping functions depends on the limiting behavior of turbu- 
lence near the wall. A correct formulation should fulfill the following relations: 

U = o(y) k = o(y 2 ) e = o(l) 

6 = o(y 2 ) n t = o(y 3 ) uv = o(y 3 ) 

together with a negative slope of the dissipation rate in the wall region. The two 

equation model proposed by Shih[23] matches the limiting forms remarkably well, 
but suffers numerical problems caused by i. In sh the decay of dissipation rate is 

-c 2 / 2 f (7) 

which is of o(l) at the wall. The turbulent viscosity is computed as 


k 2 

v t - c Re — (8) 

in which the damping function / M is based on a fourth order polynomial of y+. The 
e, needed in equations (7,8), is computed as 


? = e — 



2 k Re 


( 9 ) 


This expression gives € = o(y 2 ) at the wall. To avoid the convergence problems caused 
by unphysical overshoots of the second term on the right hand side of (9) near the 
wall (which causes negative dissipation rates), the ms model defines e as 

e = e/ t (10) 
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where the f e function depends on the turbulent Reynolds number, R t - !*Llk : 

ft = f{Rt) = 1 - exp(-y/Rt) 

At the wall, Rt = o(y 4 ) so that the function f ( is o(j/ 2 ). Equation (10) prevents 

negative dissipation rates, and makes e = o(y 2 ) at the wall. The f t function is 
designed to give e w e for y + > 6. 

Table 1 shows that all the LR forms, with the exception of jl and ms models, 
use damping functions for v t based on y or y+, while / 2 is based on R t . The use 
of y + as the exponent of the damping functions gives unphysical results in the case 
of stagnation points. In fact, the definition of the friction velocity u T appearing in 
equation (6) is: 



where n is the direction normal to the wall. At separation and reattachment points 

0 

and y+ -+ 0 regardless of the value of y. This implies that is zero all along the flow 
section making the viscous layer thickness unphysically unbounded. There are several 
known tricks to overcome this problem, like relating u r to the k peak in the viscous 
layer via the wall function for the turbulent kinetic energy, or replacing the velocity 
gradient by the maximum value of the vorticity o> in the cross flow direction. While 
these and other tricks remove the singularity, they are both arbitrary and potential 
sources of inaccuracies. Moreover, a general turbulence model should not require 
information about the flow domain geometry, like the wall distance. The authors[17] 
have shown that, although the replacement of y + is not an easy task, it is possible 
to use the ratio R L of the turbulence length scale L and the viscous length scale L„ 
defined as 



Rl approaches zero at the wall because L -> 0, L u oo, and reaches asymptotically 
a maximum in the log-law region. In the wall region the following limiting forms hold: 

Accordingly, the exponent of the damping function must be Rj to give f u = o(v) at 
the wall. The final form of is: 


/m = 1 


e^t-Ca) ' “ p H" 1 “ p ( c “ ! fli ' /4 )) 


( 12 ) 
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in which c mX and c p2 are two empirical constants. The application of ms to a fully 

developed channel flow [17] showed a much better agreement with DNS than the jl 
model. For any further detail about the models, the reader is referred to the original 
bibliography. 


3 The implicit factored solver 


The strong non-linearities usually associated with the LR forms require a robust 
implicit solver to avoid stringent stability limits. Among the several possibilities, 
the implicit approximate factorization by Beam and Warming [2] proved adequate for 
implementation in complex flows [15]. To take full advantage of the implicit time 
marching procedure for incompressible flows, the mass conservation equation, that 
does not contain time derivatives, has been modified according to the artificial com- 
pressibility approach originally developed by Chorin[5]. The corrected form of the 
continuity equation reads: 

1 dp dUi _ 

■ q: + w w — 0 

p dt ox{ 


in which p — put is an artificial equation of state and u> is the artificial compressibility 
parameter. The formulation ensures mass conservation only when a steady solution is 
reached and, consequently, ^ 0 , thereby making it impossible to follow a physical 

time transient. 

In two dimensions, the artificial compressibility, momentum, and transport equa- 
tions for e and k may be written in vector form: 


dQ dEc dFc = dE v dF v 

dt dx dy dx dy 


(13) 


where the unknown vector is Q = (p, J7, V, e, fc) T , E c , F c are convective terms, E v , F v 
axe diffusive terms, and H gathers the sink-source terms. The full forms of the 
linearization for the artificial compressibility algorithm may be found in ref. [15]. In 
order to be able to cope with any geometry, the transport equations are discretized in 
a generalized curvilinear nonorthogonal coordinate system. The differencing stencil 
is 

8 _ d_ d 

dxi ~ dC x ' + dr} Vxi 

where are the 2 — d cartesian coordinates and £, 77 are the 2 — d curvilinear coordi- 
nates. Fluxes are computed by using conservative centered finite differences. 

The implicit factored solver for a system of partial differential equations is based 
on a time linearization of the spatial operators. Since the artificial compressibility 
precludes following a physical time transient, the approximate factorization is imple- 
mented in its one step version that is first order accurate in time. The final form of 
the factored solver is 


I + 6 At 
(l + 0A t 


-atHj + - = RHS 

-a.Hj + - 0]) • A Q = A Q' 


(14) 
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in which RHS = At(-^ - ^ ^ ^ + ff), A,B are the convective Ja- 

cobians, iZ, T are the diffusive Jacobians, Hj is the source Jacobian, and A Q — 

Q t+1 — Q The solution process requires the direct inversion of block tridiagonal 
matrices. 

Particular attention was devoted to the local time step formula. A conservative 
expression was obtained by computing the eigenvalues A of the convective Jacobian: 

= U d 


^(contmutty) 


^(x- momentum) = ^ + ^ ^ 

^-momentum) = ^ ^ . Rj 

where R 2 d = d\ + dy, U d = (U d x -f V d y ) is the unsealed contravariant velocity, 
(d = £ or 77 ), and a ^ is an artificial speed of sound defined as 





The local time step formula is: 


At = 


CFL 

+ U r 7 + (R{ + Rrj ) 


(15) 


Equation (15) proved adequate for both N — S and k - e with CFL varying in 
the range 5-15. No artificial damping terms are added to the right or left hand 
side. Any stability or pressure- velocity decoupling problem was solved by refining the 
computational grid or lowering the CFL in (15). 

The time marching procedure updates the flow variables, p, {/, V by solving the 
artificial compressibility equation. The k — e equations are solved one at a time 
decoupled from the flow variables. In this manner the k — e solution lags one time 
step with respect to the N — S equations. 


4 Linearization of source terms 


The Hj Jacobian of the source term vector iT, that 
as 


Hj = 


oh 

dQ 


appears in equation (14) is defined 


While the derivation of the convective and diffusive Jacobians is a straightforward 
exercise, the choice of the linearization of H may play a significant role in the ro- 
bustness of a turbulent flow solver. H can not be treated in a fully explicit way and, 
to maintain the stability of iterative procedure, it must be properly linearized. Hj is 
usually introduced in the sweep associated with the largest gradients, but in elliptic 
flows this direction is not always unique, so that it was found convenient to introduce 
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model 

dH s 

de 

dHj. 

~dk^ ..... 

ch 

- 2 / 2 I-! 

2 

Re y 2 

ji 

-C 2 h\ 

e 

k 

nh 

-C2 

€ 

k 

CO 

- 2c 4~2 

€ 

2 k 

sp 

■-s* p 

€ 

k 

ki 

- 2 c 2 £ ( for both e eq.) 

— | ( for both k eq.) 

TO 

-2C2^ 

€ 

! k 

lb 

-2<^ 

£ 

! k 

sh 

" ‘2/2f 

1 € 
k 

ms 

. “ C2 M 

\ € 

L b 


Table 3: Approximate source Jacobian 


it in both the sweeps according to the directional weights a in equation (14). Typical 
values of a are 1 — 5. 

The source Jacobian i/j, if negative, improves the algorithm robustness by in- 
creasing the diagonal dominance of the tridiagonal block matrix. This observation 
indicates that the optimal form of Hj is not necessarily the exact differentiation of 
H , that may have positive and negative contributions, but an approximate form that 
is always negative. Possible numerical problems may be also caused by overshoots of 
Hj. In fact the tridiagonal matrix given by equation (14) is 

M-i • A + (Hj + A{) • AUi + A i+l • A U i+l - RHS { 

If H j > A and Hj > RHS , A U — ► 0 regardless of the value of RHS. This may 
unphysically freeze the local solution. 

After intense numerical testing, we found it convenient to compute Hj dropping 
the positive terms. To prevent possible over- or under-shoots only the leading terms 
in the viscous layer were retained. They have been determined by substituting in 
fc, e, U their limiting forms at the solid boundary for each model. The linearization 
for the ten models may be found in table 3. The analysis of the source Jacobians of 
the formulations showed that all the models used in this paper have the same degree 
of the robustness with the exception represented by the lb model for which serious 
numerical problems were encountered. 


5 Turbulent flow past a hill 

The aforedescribed models have been compared for a fully developed channel flow 
[17], The computed results showed significant discrepancies among the models and 
only few of the formulations proved to fit with the DNS data at Re = 3300. The very 
low Reynolds number and the simple geometry used for testing indicated that further 
investigation is necessary to verify the model performances in a complex geometry at 
very high iZe, as often encountered in engineering applications. 
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grid # 

points 

st 

1 

111 x 81 

1.08 

2 

125 x 121 

1.07 

3 

221 X 121 

1.07 


Table 4: Grids (st^stretch ratio) 


For incompressible two-dimensional flows with streamline curvature there are only 
a few valuable experimental data sets. Among these we chose the flow past a circular 
hill[l]. In this test an incoming boundary layer experiences a short region of concave 
curvature with adverse pressure gradient, prior to encountering a long region of con- 
vex curvature in which the pressure gradient reverts to favorable and then becomes 
adverse. At this point the measurements detect a small separated region downstream 
of the hill in correspondence to a short convex curvature of the lower boundary. This 
experiment provides an interesting opportunity to test the turbulence models in a 
mixed adverse-favorable pressure gradient condition. The flow approaches the hill 
with a thin boundary layer (£ = 55.7 mm). The reference velocity U re f is approx- 
imately 20 m/s, which gives a Reynolds number based on the hill curvature radius 
(= 1284 mm) approximately equal to 1.4 • 10 6 . 

The measurements include the skin friction coefficient, c/, obtained by the Clauser 
chart and the Preston tube all along the hill together with the wall static pressure 
distribution, c p . The two methods used to determine Cf did not give significantly 
different results. 

The mean velocity, turbulent normal and shear stresses are measured in ten dif- 
ferent sections along the hill. The mean velocities are measured by a Pitot tube 
accounting for displacement effects in the wall proximity. The turbulent stresses are 
measured by using normal and cross hot-wires. The cross sectional measurements are 
interrupted very close to separation: this prevents a detailed comparison in the region 
well inside the recirculation. 

The computations have been carried out using the three grids summarized in table 
4; figure 1 shows grid number 2. The stretch ratio st is defined as the ratio of two 
consecutive grid cell heights in the cross flow direction. 

The grid height was twice the hill curvature radius to minimize the effect of 
the upper boundary on the solution, where a zero normal gradient condition was 
implemented. On the exit boundary, a zero gradient condition was found adequate 
provided that the computational domain was extended sufficiently downstream of the 
hill. 

Because of the incoming boundary layer and the large Reynolds number, grid 
number 1 allowed placing the first grid point away from the wail at y + zz 9, while 
with grids number 2 and 3 the points were more clustered at the wall to give % 1. 
Grid 3 was introduced to check whether a stronger mesh refinement could bring about 
further changes in the computed results. The computations showed that the converged 
results obtained with grid 2 were nearly indistinguishable from those obtained with 
grid 3, so that all the final calculations were carried out with grid 2. 

The mean velocity and turbulence quantities at the inlet section of the compu- 
tational domain were specified according to the experimental profiles at s = 596. 
An attempt to use a fully developed flow inlet condition revealed the extreme sen- 
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sitivity of the computed results to the incoming profiles. When the inlet boundary 
layer thickness or the turbulent kinetic energy did not match with measurements, the 
computed profiles were far off from experiments. Figure 2 shows the computed tur- 
bulent viscosity levels at X = 1990 mm for both the fully developed and experimental 
inlet conditions. The v t profiles differ by approximately one order of magnitude im- 
mediately outside of the buffer layer. This observation shows the stringent need for 
detailed inlet conditions to perform a meaningful computation. 


5,1 Wall static pressure 

First let us look at the static pressure, <7 P , along the lower side of the computational 
domain. Figure 3 gives the computed C p using grids 1 and 2 with various models. 

The ch model shows a good fit with experiments in the adverse pressure gradient 
(apg) region corresponding to the positive peak in C p . The peak is well captured by 
both the grids. The same degree of accuracy is maintained all along the favourable 
pressure gradient ( fpg ) region. When the gradient reverts to positive the two grids 
start to give different results. The coarse grid generally produces larger recirculations 
with a moderate pressure recovery downstream of the hill, but it completely fails 
to reproduce the other flow variables, whereas the refined grid number 2 predicts a 
stronger pressure recovery than that shown by experiments. Nearly all the models 
gave a shorter and thinner recirculation with the refined grid as compared to the 
coarse one. This is somewhat surprising, but the converged results clearly showed 
that a large reversed flow region was detected by using the coarse grid number 1 
together with very thick boundary layers in contrast with experiments in which the 
boundary layer remains quite thin. The same excessive pressure recovery was also 
found by Kim[8]. 

The lb model seems to follow the experimental C p distribution quite closely and 
clearly predicts the static pressure drop in the large separated region downstream of 
the hill. The grid independence is quite remarkable. Unfortunately this model caused 
serious numerical problems and the results with both grids 1 and 2 are only partially 
converged. 

The nh model is very similar to the jl formulation. However, the two formulations 
give results considerably different for a fully developed channel flow at moderate Re. 
Figure 3 shows that the differences between the models fade out for a high Re flow 
provided that a refined grid is used. In fact when using a coarse grid the inaccuracies 
in the k and e profiles directly affect the turbulent viscosity distribution via the 
function based on Rt. The same effect is clearly smeared when using based on 
which is much less sensitive to the local turbulent flow field. The differences between 
the coarse and the refined grid are larger in the adverse pressure gradient region in 
which the flow is likely to encounter the largest gradients. The direct comparison 
of these two models shows that, even if the replacement of the damping functions 
based on is desirable, the use of Rt introduces a strong mesh dependence and the 
use of a very refined grid becomes compulsory to obtain a reasonable static pressure 
distribution (typically, the first grid point away from the wall must be placed at 
t/+ ~ 1, as in grid 2, to have grid independent results). 

The q — u; model introduced by Coakley, co, was found to give a poor fit with DNS 
of a channel flow[17]. Nevertheless, it succeeds in giving a satisfactory agreement with 
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the measured flow past a hill. The model appears to be very grid sensitive on account 
of the differences in the static pressure distribution in both the apg regions. The 
main reason for this is again the damping function which is based on R y . As in 
the jl model, if the peaks of k and e are not properly resolved, the inaccuracies will 
be felt also in the damping of the turbulent viscosity with large effects on the mean 
flow field. 

The sp model fits remarkably well the experimental static pressure coefficient 
until the separated region is approached. Moreover the model seems relatively grid 
independent until the strong apg is reached. At this point sp behaves like the other 
formulations. The pressure recovery in the separated region is still too large, and, 
when comparing this model with the other LR forms that use the wall distance t/, it 
is evident that sp does not bring about any considerable improvement. 

The sh and ms forms are very similar and, when using the refined grid, the 
static pressure profiles are nearly indistinguishable. The introduction of the pressure 
transport term II does not seem to play a significant role in the determination of the 
static pressure. For these two models the effect of the damping function / M based on 
y+ or R t is even more evident. Let’s recall that in sh the function is based on 
a fourth order polynomial of y+, while ms uses a complex function of the turbulent 
Reynolds number. Since / M appears in both the turbulent viscosity formula and the 
pressure transport term, the ms model will be very grid dependent for the same reason 
that jl was found to be grid dependent. This feature of ms was not highlighted in 
the channel flow computations [17] because of the much less stringent mesh points 
requirement. 

The two-layer models give results that are somewhat less grid sensitive than the 
other formulations. This is partly due to the algebraic dissipation rate equation in the 
wall proximity, which, if solved by a transport equation down to the wall, appears to 
have more stringent mesh requirements. Nevertheless ro and ki predictions are very 
similar to the other models in the first apg and fpg regions. Again, differences start 
to appear in the strong apg region prior to separation. Here, while ro shows a notable 
difference between the two grids, ki gives profiles remarkably grid independent. It 
seems that the two scale model serves to relax the grid requirement. 


5.2 Skin friction 

The skin friction coefficient, nondimensionalized with respect to the free stream dy- 
namic pressure, exhibits a peak approximately corresponding to the top of the hill. 
It reverts to negative values downstream of the hill in the short region with con- 
vex curvature where separation takes place. The plots report only the skin friction 
coefficient, Cf , computed with the refined grid number 2. 

The form of the friction velocity u r to be used in the damping functions based 
on j/+ is of primary importance for the skin friction coefficient. To avoid numerical 
problems caused by sudden changes in the turbulent viscosity, the friction velocity 
(mentioned in section 2.1) was related to the peak of turbulent kinetic energy closest 
to the wall via the standard wall function for k 
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Small convex curvature can produce a significant reduction of the skin friction 
coefficient, whereas a concave curvature causes an increase in Cf. In the present 
case, the ratio 8/R, where R is the local curvature radius, is above 0.01 so that the 
curvature effect gives strong changes in the skin friction. 

Figure 5 reports the comparison between the measured and computed distribu- 
tion. The ch formulation shows a satisfactory agreement with the experimented data 
in both the apg and fpg regions. In the first apg region ch seems to reproduce the 
measurements quite well. The sudden inversion from apg to fpg gives a very thin 
boundary layer together with a dramatic increase in the skin friction. The behavior 
of Cf quite closely resembles that of the static pressure coefficient. The peak in Cf 
is located approximately at the static pressure minimum prior to inversion to apg . 
The separation point is not fixed by geometrical constraints and it results from the 
local boundary layer thickness and pressure gradient. The general tendency of the 
models was to locate the separation point downstream of the measured one which 

is at X[ exp) % 2095 mm. ch locates the separation at X B % 2155 mm. The profile 
given by lb is fax from the experimental results. Moreover lb locates a short recircu- 
lation upstream of the hill corresponding to a strong decelerated region at the static 
pressure maximum in contrast with measurements where no separation is detected. 
Nevertheless, the lb results, because of the aforementioned numerical problem, are 
only partially converged. 

The jl and nh models give similar results. The separation point is located by 
both the models at X B % 2170 mm. The direct comparison of the two models allows 
evaluating the effect of the damping function based on R t or on y+. The low Re 
calculations of ref.[17] showed that the damping form selected in jl was suited only 
for high Re flows. In fact, in this test, it appears that the boundary layer is properly 
predicted together with the shape of the skin friction coefficient. On the other hand, 
the damping function in nh which requires does not perform as well as the one 
in especially in the fpg region. 

The formulation proposed by Coakley (co) is the only one that does not detect any 
reversed flow downstream of the hill. The DNS comparison for the channel flow[17] 
showed that co clearly overestimates the turbulent viscosity all along the cross section. 
This ends in producing excessive momentum diffusion and less or no separation at all. 
The sp model follows the measurement remarkably well in the fpg region and locates 
the separation point at approximately X B % 2150 mm. sp responds to apg with a 
dramatic drop in the skin friction. While this behavior is qualitatively correct, the 
Cf drop in both the apg regions is steeper than the experimental drop. This feature 
is in contrast with all the other models in which there is a tendency to closely follow 
the experimental skin friction. 

The two formulations that account for the pressure transport term, II, behave 
differently on account of the selected form of the damping function f^. While sh gives 
probably the best agreement with experiments, ms overestimates the skin friction 
peak by approximately 5%. The pressure transport term is more important in apg 
and, in fact, sh and ms seem to give an improved fit with experiments as soon as the 
pressure gradient reverts from favorable to adverse. The introduction of the pressure 
correlation in the transport of the turbulent kinetic energy proves here to have a 
beneficial effect on the prediction. Although further testing is necessary to draw a 
final conclusion, it seems that the modeling of II as a turbulent kinetic energy diffusion 
is appropriate at least for this class of flows. 
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Table 5: Exponents of the wall limiting forms* 


For the two-layer models the effect of the matching criteria on Cf was investigated. 
Figure 6 shows the interface distance from the wall for both to and ki. The two 
models predict different turbulent viscosity levels in the buffer layer. In fact, while the 
automatic matching criteria was i/t ~ 36 for both models, this value of the turbulent 
viscosity is reached at a larger distance from the wall by ki than to . Nevertheless, 
the difference is at most 0.1% of the hill curvature radius. Moreover this difference 
seems to be constant all along the hill shape and it qualitatively resembles the shape 
of the wall static pressure for both the models. The to formulation fails to predict 
separation when fixing the matching point in the computational domain according to 
the experimental inlet profiles. Conversely, the separation point is remarkably well 
predicted at X, % 2095 mm when letting the model adjust the matching level at v t « 
36 in accordance to the locally computed turbulent viscosity. The to formulation has 
already proved to give larger separation bubbles with respect to other LR forms[14], 
and this test confirms this tendency. The reason for this is probably that in ro the 
dissipation rate is computed by using the Norris and Reynolds algebraic expression 
within the buffer layer which gave better results than the standard k - 6 model in 
apg[ 22], The results given by ki are similar. 


5.3 Cross sections 

In order to provide a complete view of the model performance, the mean velocity, tur- 
bulent kinetic energy, and turbulent shear stress profiles are reported in five different 
cross sections. The location of the five sections may be found in Fig.l. The results 
reported in this set of plots refer to grid number 2 only. 

In table 5 all the ten models show similar limiting forms at the wall for the mean 
velocity and turbulent kinetic energy. Conversely, the turbulent shear stress modeled 
by nh, co, and lb is o(y 4 ) at the wall, whereas the other models show o(y 3 ). 

Figure 7 reports the mean velocity profiles at the aforementioned five cross sec- 
tions, while figures 8 and 9 show the turbulent shear stress and the turbulent kinetic 
energy respectively. The plots show a closeup view of the region close to the wall 
since the flow variables were nearly flat in the outer part of the domain. 

The first section at x = 596 mm is not of particular interest since the disturbance 
given by the presence of the hill is not yet felt. Nevertheless, all the models agree with 
experiments within a narrow range, especially for velocity and turbulent shear stress 
profiles. It is observed that there is no tZ^tZj inversion and that the turbulent shear 
stress reaches zero at Y « 60 for all the models. The computed turbulent kinetic 
energy is underestimated by all the models. 

The second section at X = 1183 mm is of much more interest since it is located 
near the pressure gradient inversion immediately downstream of the beginning of 
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the long convex region. Both ch and lb fit the measured velocity profiles which 
show a very thin boundary layer caused by the presence of the obstacle, jl and nh 
seem to overestimate the boundary layer thickness, (jl produces a much stronger 
overestimation of the boundary layer thickness than nh). 

The strong reduction of the boundary layer thickness results in a large mean strain 
in the viscous layer. This can be seen in the peak of the turbulent shear stress and 
turbulent kinetic energy located in the wall proximity, jl and nh show a strong k 
overshoot that is not found by the other models. In checking every single model term, 
this behavior appeared to be caused by the selected form of the extra destruction rate, 
l?, that is proportional to the square of the turbulent kinetic energy gradient as follows 


D cc 



(16) 


Equation (16) seems to produce k overshoots in fpg regions. 

The third cross flow section, at X — 1469 mm is well inside the fpg region. With 
respect to the previous section the growth of the boundary layer is evident even if 
the overall velocity profile seems to be only weakly affected. The reduction of the 
turbulent shear stress is well predicted by all the models, jl and nh give overshoots 
in k and u^u], and qualitatively behave as in the previous section. 

When moving to the fourth and fifth sections at X = 1862 mm and X - 1990 mm 
respectively, differences between the model predictions start to be more evident. In 
these sections the flow encounters a much stronger apg than that found before section 
2. The predicted velocity profiles are in good agreement with measurements for 
nearly all the models. Nevertheless the velocity profiles predicted by sh and ms seem 
to be in slightly better agreement with experiments than those predicted by using the 
other models, especially at the edge of the boundary layer. All the models tend to 
overestimate ujuj although the shape of the measured turbulent shear stress is well 
predicted. 


5.4 Numerical details 

Sections 3 and 4 give some details about the solver structure that was left unchanged 
for all the models. In the calculations the turbulence field was implicitly updated at 
the end of every implicit N — S sweep. No advantages were found in updating k and 
e after a fixed number of N — S sweeps. The computations were carried out on the 
Cray YMP computer at NASA Lewis Research Center. Converged results with the 
125 X 121 grid were obtained after 2000 iterations in % 500 - 600 seconds. The slow 
convergence rate was caused by the lack of numerical damping terms and by the grid 
cell aspect ratio that reached several thousands. 

The only LR form showing serious numerical instabilities was lb. The proposed 
linearization for the source terms, while preventing the code from blowing up, did not 
achieve the overall residual of the order of 10~ 6 - 10 -7 that was requested for all the 
other models. 
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6 Concluding Remarks 

The comparison of ten different LR two-equation models in a complex geometry 
shows the limits of these models. Despite the improvements in channel flow pre- 
dictions (highlighted in a previous investigation[17]) introduced by the use of the 
low-Reynolds-number corrections to the high- Reynolds number two-equation model, 
problems occur when trying to model flows with strong pressure gradients and stream- 
line curvature. The large differences among the model predictions for a simple channel 
flow at a moderate Reynolds number seem to fade out in such a complex geometry. 
The reason for that is not very clear. This may be, in part, due to the very high 
Reynolds number that reduces the overall influence of the LR model corrections, sh , 
and consequently ms, attempts to introduce an extra model term which is the pres- 
sure correlation II in the k equation. The skin friction coefficient seems to be the 
only quantity affected by II in response to the pressure gradients and the other flow 
variables appear to be only marginally influenced. 

In addition, the ms model does not explicitly depend on the wall distance and 
yields comparable results with the best of the other two-equation models. This indi- 
cates its promise for modeling of complex turbulent flows. 

When using a single solver to carry out calculations with different LR models the 
differences in the predictions must be attributed only to the model formulations. Con- 
sequently, the present set of tests gives some understanding of the model discrepancies 
that were hidden behind purely numerical differences. 
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Figure 3: Wall static pressure 
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Figure 5: Skin friction coefficient (contd) 
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Figure 6: Match level for two-layer models 


25 






























X = 1469 X- 1469 


I 

■4 


X= 1469 


0.0 -c 

).5 

G.O 

-0.5 

0.0 

-0.5 

0.0 

-0.5 

uv* tOO 


uvMOO 


uv* 100 


o 

o 

uv* 100 

X=1862 


X= 1862 


X= 1862 


X=1862 

X=1862 


uv*100 


Figure 8: Cross sections: turbulent shear stress (contd) 
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